using System;
using dnAnalytics.Nli;

/// <summary>
/// LU decompisition.
/// </summary>
internal class Program {
	[STAThread]
	private static void Main() {
		//use the MKL LAPACK provider
		ILapack lapack = Lapack.Provider(LapackProvider.Mkl);
		IBlas blas = Blas.Provider(BlasProvider.Mkl);

		/*
		    1 2 3             -0.625  0.458 -0.167
		a = 5 6 7   inv(a) =  -0.500 -0.167  0.333
			4 9 8              0.875 -0.042 -0.167
		*/

		//row ordering
		double[] a = new double[] {1, 5, 4, 2, 6, 9, 3, 7, 8};
		double[] inv = new double[] {1, 5, 4, 2, 6, 9, 3, 7, 8};

		//note: the indices ipiv contains are 1 based, not 0 based.
		int[] ipiv;
		lapack.Dgetrf(3, 3, inv, 3, out ipiv);
		lapack.Dgetri(3, inv, 3, ipiv);

		foreach (double value in inv) {
			Console.WriteLine(value.ToString("0.###"));
		}
		Console.WriteLine();

		//verify the a * inv(a) = identity.
		double[] ident = new double[] {0, 0, 0, 0, 0, 0, 0, 0, 0};
		blas.Dgemm(BlasOrderType.Column, BlasTransType.NoTrans, BlasTransType.NoTrans, 3, 3, 3, 1, a, 3, inv, 3, 0, ident, 3);
		foreach (double value in ident) {
			Console.WriteLine(value.ToString("0.###"));
		}		
		
		/* Output:
		-0.625
		-0.5
		0.875
		0.458
		-0.167
		-0.042
		-0.167
		0.333
		-0.167

		1
		0
		0
		0
		1
		0
		0
		0
		1
		*/
	}
}